#include "context.h"
#include "evaluator_utils.h"
#include <stdint.h>
using namespace NTL;
uint64_t Context::BASECONVERT = 0;
/*
 * @brief Construct a new Context:: Context object
 *
 * @param logN:多项式次数
 * @param L: 模数Q长度
 * @param K: 模数P长度
 * @param T: 模数T长度
 * @param h: 秘钥稀疏度
 * @param sigma: 高斯分布标准差
 */
Context::Context(long logN, long L, long K, long T, int lambda, int dNum, long h, double sigma)
    : logN(logN), L(L), K(K), T(T), h(h), sigma(sigma)
{
    if(lambda == 0) {
        if ((logN == 15) && (L == 15)) {
            GetParameters(params::logN15QP881);
        }
        else if ((logN == 12) && (L == 12)) {
            GetParameters(params::logN12QP884);
        }
        else {
            cout << "没有满足要求的参数！！" << endl;
        }
        if (L > std::numeric_limits<size_t>::max() / sizeof(uint64_t)) {
            return;
        }
        if (K > std::numeric_limits<size_t>::max() / sizeof(uint64_t)) {
            return;
        }
        if (T > std::numeric_limits<size_t>::max() / sizeof(uint64_t)) {
            return;
        }
    }
    else if ((lambda == 128) || (lambda == 192)){
        logp = 55;
        GetParametersLambda(logN, lambda, dNum);
    }
    else {
        throw std::domain_error("Invalid value for lambda.");
    }

    L = this->L;
    K = this->K;
    N = 1L << logN;
    M = N << 1;
    logNh = logN - 1;
    Nh = N >> 1;
    BASECONVERT = L; 
    ex = new uint64_t[(L + K) << logN];
    vx = new uint64_t[(L + K) << logN];
    axt = new uint64_t[L << logN];
    bxt = new uint64_t[L << logN];
    axbx1 = new uint64_t[L << logN];
    axbx2 = new uint64_t[L << logN];
    axax = new uint64_t[L << logN];
    bxbx = new uint64_t[L << logN];
    tmpAx = new uint64_t[(L + K) << logN];
    tmpBx = new uint64_t[(L + K) << logN];
    cwQ = new uint64_t[L << logN];
    cwP = new uint64_t[K << logN];
    new uint64_t[K << logN];
    cwNTT = new uint64_t[L << logN];
    cwQSet = new uint64_t*[L];
    cwPSet = new uint64_t*[L];
    #pragma omp parallel for
    for(int i = 0;i < L; i++) {
        cwQSet[i] = new uint64_t[L * N]();
        cwPSet[i] = new uint64_t[K * N]();
    }
    polAx = new uint64_t[(L + K) << logN];
    polBx = new uint64_t[(L + K) << logN];
    cxInvNTT = new uint64_t[L << logN];
    tmpCtAx = new uint64_t[L << logN];
    tmpCtBx = new uint64_t[L << logN];
    tmpP2Ql = new uint64_t[K << logN];
    c2QP = new uint64_t[(L + K) << logN];
    tmpQl2QlP = new uint64_t[L << logN];
    b = new uint64_t[N];
    tmpQ2T = new uint64_t[L << logN];
    twoqVec = new uint64_t[L]();
    qrVec = new __uint128_t[L];
    qrVec_high = new uint64_t[L]();
    qrVec_low = new uint64_t[L]();
    qdVec = new uint64_t[L];
    qInvVec = new uint64_t[L];
    qRoots = new uint64_t[L];
    qRootsInv = new uint64_t[L];
    qRootPows = new uint64_t *[L];
    qRootPowsInv = new uint64_t *[L];
    NInvModq = new uint64_t[L];
    qRootPows_shoup = new uint64_t *[L];
    qRootPowsInv_shoup = new uint64_t *[L];
    NInvModq_shoup = new uint64_t[L];
    for (long i = 0; i < L; i++) {
        qrVec[i] = ~static_cast<__uint128_t>(0) / qVec[i];
        twoqVec[i] = static_cast<uint64_t>(2 * qVec[i]);
        qrVec_high[i] = static_cast<uint64_t>(qrVec[i] >> 64);
        qrVec_low[i] = static_cast<uint64_t>(qrVec[i]);
        qdVec[i] = qVec[i] << 1;
        qRoots[i] = FindMthRootOfUnity(M, qVec[i]);
        qRootsInv[i] = InvMod(qRoots[i], qVec[i]);
        NInvModq[i] = InvMod(static_cast<uint64_t>(N), qVec[i]);
        qInvVec[i] = Inv(qVec[i]);
        qRootPows[i] = new uint64_t[N];
        qRootPowsInv[i] = new uint64_t[N];
        uint64_t power = static_cast<uint64_t>(1);
        uint64_t powerInv = static_cast<uint64_t>(1);
        qRootPows_shoup[i] = new uint64_t[N]();
        qRootPowsInv_shoup[i] = new uint64_t[N]();
        NInvModq_shoup[i] = x_Shoup(NInvModq[i], qVec[i]);
        for (long j = 0; j < N; j++) {
            uint64_t jPrime = BitReverse(static_cast<uint32_t>(j)) >> (32 - logN);
            qRootPows[i][jPrime] = power;
            qRootPows_shoup[i][jPrime] = x_Shoup(power, qVec[i]);
            unsigned __int128 tmp = (static_cast<unsigned __int128>(power) << 64);
            qRootPowsInv[i][jPrime] = powerInv;
            qRootPowsInv_shoup[i][jPrime] = x_Shoup(powerInv, qVec[i]);
            if (j < N - 1) {
                MulMod(power, power, qRoots[i], qVec[i]);
                MulMod(powerInv, powerInv, qRootsInv[i], qVec[i]);
            }
        }
    }
    twopVec = new uint64_t[L]();
    prVec_high = new uint64_t[L]();
    prVec_low = new uint64_t[L]();
    prVec = new __uint128_t[L];
    pdVec = new uint64_t[K];
    pInvVec = new uint64_t[K];
    pRoots = new uint64_t[K];
    pRootsInv = new uint64_t[K];
    pRootPows = new uint64_t *[K];
    pRootPowsInv = new uint64_t *[K];
    NInvModp = new uint64_t[K];
    NInvModp_shoup = new uint64_t[K];
    pRootPows_shoup = new uint64_t *[K];
    pRootPowsInv_shoup = new uint64_t *[K];
    for (long i = 0; i < K; i++) {
        prVec[i] = ~static_cast<__uint128_t>(0) / pVec[i];
        twopVec[i] = static_cast<uint64_t>(2 * pVec[i]);
        prVec_high[i] = static_cast<uint64_t>(prVec[i] >> 64);
        prVec_low[i] = static_cast<uint64_t>(prVec[i]);
        pdVec[i] = pVec[i] << 1;
        pRoots[i] = FindMthRootOfUnity(M, pVec[i]);
        pRootsInv[i] = InvMod(pRoots[i], pVec[i]);
        NInvModp[i] = InvMod(static_cast<uint64_t>(N), pVec[i]);
        pRootPows[i] = new uint64_t[N];
        pRootPowsInv[i] = new uint64_t[N];
        pInvVec[i] = Inv(pVec[i]);
        uint64_t power = static_cast<uint64_t>(1);
        uint64_t powerInv = static_cast<uint64_t>(1);
        pRootPows_shoup[i] = new uint64_t[N]();
        pRootPowsInv_shoup[i] = new uint64_t[N]();
        NInvModp_shoup[i] = x_Shoup(NInvModp[i], pVec[i]);
        for (long j = 0; j < N; j++) {
            uint64_t jPrime = BitReverse(static_cast<uint32_t>(j)) >> (32 - logN);
            pRootPows[i][jPrime] = power;
            pRootPows_shoup[i][jPrime] = x_Shoup(power, pVec[i]);
            unsigned __int128 tmp = (static_cast<unsigned __int128>(power) << 64);
            pRootPowsInv[i][jPrime] = powerInv;
            pRootPowsInv_shoup[i][jPrime] = x_Shoup(powerInv, pVec[i]);
            if (j < N - 1) {
                MulMod(power, power, pRoots[i], pVec[i]);
                MulMod(powerInv, powerInv, pRootsInv[i], pVec[i]);
            }
        }
    }
    qHatModq = new uint64_t *[L];
    for (long l = 0; l < L; l++) {
        qHatModq[l] = new uint64_t[l + 1];
        for (long i = 0; i < l + 1; i++) {
            qHatModq[l][i] = 1;
            for (long j = 0; j < i; j++) {
                uint64_t temp = qVec[j] % qVec[i];
                MulMod(qHatModq[l][i], qHatModq[l][i], temp, qVec[i]);
            }
            for (long j = i + 1; j < l + 1; j++) {
                uint64_t temp = qVec[j] % qVec[i];
                MulMod(qHatModq[l][i], qHatModq[l][i], temp, qVec[i]);
            }
        }
    }
    qHatInvModq = new uint64_t *[L];
    for (long l = 0; l < L; l++) {
        qHatInvModq[l] = new uint64_t[l + 1];
        for (long i = 0; i < l + 1; i++) {
            qHatInvModq[l][i] = InvMod(qHatModq[l][i], qVec[i]);
        }
    }
    PModq = new uint64_t[L];
    for (long i = 0; i < L; i++) {
        PModq[i] = 1;
        for (long k = 0; k < K; k++) {
            uint64_t temp = pVec[k] % qVec[i];
            MulMod(PModq[i], PModq[i], temp, qVec[i]);
        }
    }
    QModp = new uint64_t *[L];
    for (long i = 0; i < L; i++) {
        QModp[i] = new uint64_t[K];
        for (long k = 0; k < K; k++) {
            QModp[i][k] = 1;
            for (long j = 0; j < i + 1; j++) {
                uint64_t temp = qVec[j] % pVec[k];
                MulMod(QModp[i][k], QModp[i][k], temp, pVec[k]);
            }
        }
    }
    qInvModq = new uint64_t *[L];
    for (long i = 0; i < L; i++) {
        qInvModq[i] = new uint64_t[L];
        for (long j = 0; j < i; j++) {
            qInvModq[i][j] = InvMod(qVec[i], qVec[j]);
        }
        for (long j = i + 1; j < L; j++) {
            qInvModq[i][j] = InvMod(qVec[i], qVec[j]);
        }
    }
    rotGroup = new long[Nh];
    long fivePows = 1;
    for (long i = 0; i < Nh; i++) {
        rotGroup[i] = fivePows;
        fivePows *= 5;
        fivePows %= M;
    }
    trVec = new __uint128_t[T];
    twotVec = new uint64_t[L]();
    trVec_high = new uint64_t[L]();
    trVec_low = new uint64_t[L]();
    tInvVec = new uint64_t[T];
    tdVec = new uint64_t[T];
    tInvVec = new uint64_t[T];
    tRoots = new uint64_t[T];
    tRootsInv = new uint64_t[T];
    tRootPows = new uint64_t *[T];
    tRootPowsInv = new uint64_t *[T];
    NInvModt = new uint64_t[T];
    tRootPows_shoup = new uint64_t *[T];
    tRootPowsInv_shoup = new uint64_t *[T];
    NInvModt_shoup = new uint64_t[T];
    for (long i = 0; i < T; i++) {
        trVec[i] = ~static_cast<__uint128_t>(0) / tVec[i];
        twotVec[i] = static_cast<uint64_t>(2 * tVec[i]);
        trVec_high[i] = static_cast<uint64_t>(trVec[i] >> 64);
        trVec_low[i] = static_cast<uint64_t>(trVec[i]);
        tdVec[i] = tVec[i] << 1;
        tRoots[i] = FindMthRootOfUnity(M, tVec[i]);
        tRootsInv[i] = InvMod(tRoots[i], tVec[i]);
        NInvModt[i] = InvMod(static_cast<uint64_t>(N), tVec[i]);
        tRootPows[i] = new uint64_t[N];
        tRootPowsInv[i] = new uint64_t[N];
        tInvVec[i] = Inv(tVec[i]);
        uint64_t power = static_cast<uint64_t>(1);
        uint64_t powerInv = static_cast<uint64_t>(1);
        tRootPows_shoup[i] = new uint64_t[N];
        tRootPowsInv_shoup[i] = new uint64_t[N];
        NInvModt_shoup[i] = x_Shoup(NInvModt[i], tVec[i]);
        for (long j = 0; j < N; j++) {
            uint64_t jprime = BitReverse(static_cast<uint32_t>(j)) >> (32 - logN);
            tRootPows[i][jprime] = power;
            tRootPows_shoup[i][jprime] = x_Shoup(power, tVec[i]);
            unsigned __int128 tmp = (static_cast<unsigned __int128>(power) << 64);
            tRootPowsInv[i][jprime] = powerInv;
            tRootPowsInv_shoup[i][jprime] = x_Shoup(powerInv, tVec[i]);
            if (j < N - 1) {
                MulMod(power, power, tRoots[i], tVec[i]);
                MulMod(powerInv, powerInv, tRootsInv[i], tVec[i]);
            }
        }
    }
    uint64_t pow = 1, pos = 0, mask = 2 * N - 1;
    indexMatrix = new uint64_t[N];
    for (int i = 0, j = int(N >> 1); i < int(N >> 1); i++, j++) {
        pos = BitReverse64(pow >> 1, logN);
        indexMatrix[i] = pos;
        indexMatrix[j] = N - pos - 1;
        pow *= GaloisGen;
        pow &= mask;
    }
    qProduct = new ZZ[L];
    qProduct[0] = conv<ZZ>(qVec[0]);
    for (int i = 1; i < L; i++) {
        qProduct[i] = qProduct[i - 1] * conv<ZZ>(qVec[i]);
    }
    t0InvModQ = new ZZ[L];
    for (int i = 0; i < L; i++) {
        t0InvModQ[i] = InvMod(conv<ZZ>(tVec[0]), qProduct[i]);
    }
    pModulusAtLevel = new ZZ[K];
    pModulusAtLevel[0] = conv<ZZ>(pVec[0]);
    for (int i = 1; i < K; i++) {
        pModulusAtLevel[i] = pModulusAtLevel[i - 1] * conv<ZZ>(pVec[i]);
    }
    qHatModt = new uint64_t **[L];
    for (long l = 0; l < L; l++) {
        qHatModt[l] = new uint64_t *[l + 1];
        for (long i = 0; i < l + 1; i++) {
            qHatModt[l][i] = new uint64_t[T];
            for (long k = 0; k < T; k++) {
                qHatModt[l][i][k] = 1;
                for (long j = 0; j < i; j++) {
                    uint64_t temp = qVec[j] % tVec[k];
                    MulMod(qHatModt[l][i][k], qHatModt[l][i][k], temp, tVec[k]);
                }
                for (long j = i + 1; j < l + 1; j++) {
                    uint64_t temp = qVec[j] % tVec[k];
                    MulMod(qHatModt[l][i][k], qHatModt[l][i][k], temp, tVec[k]);
                }
            }
        }
    }
    QModt = new uint64_t *[L];
    for (long i = 0; i < L; i++) {
        QModt[i] = new uint64_t[T];
        for (long k = 0; k < T; k++) {
            QModt[i][k] = 1;
            for (long j = 0; j < i + 1; j++) {
                uint64_t temp = qVec[j] % tVec[k];
                MulMod(QModt[i][k], QModt[i][k], temp, tVec[k]);
            }
        }
    }
    qHatInvModqShoup = new uint64_t *[L];
    for (auto j = 0; j < L; j++) {
        qHatInvModqShoup[j] = new uint64_t[j + 1];
        for (auto k = 0; k < j + 1; k++) {
            qHatInvModqShoup[j][k] = x_Shoup(qHatInvModq[j][k], qVec[k]);
        }
    }
    pow2Base = 0;
    matrixValues = new uint64_t *[L];
    for (int i = 0; i < L; i++) {
        matrixValues[i] = new uint64_t[3];
        matrixValues[i][0] = 0;
        matrixValues[i][1] = 1;
        matrixValues[i][2] = qVec[i] - 1;
    }
    qInvModqWithQP = new uint64_t ***[K - 1];
    qInvModqpWithQP = new uint64_t ****[K - 1];
    QModqpWithQP = new uint64_t ***[K - 1];
    vector<int> xnbPi;
    vector<uint64_t> QiWithQP;
    vector<uint64_t> PiWithQP;
    for (int lvlP = 0; lvlP < K - 1; lvlP++) {
        int nbPi = lvlP + 2;
        int decompRNS = int(ceil((L * 1.0) / nbPi));
        xnbPi.resize(decompRNS);
        for (int i = 0; i < decompRNS; i++) {
            xnbPi[i] = nbPi;
        }
        if (L % nbPi != 0) {
            xnbPi[decompRNS - 1] = L % nbPi;
        }
        qInvModqWithQP[lvlP] = new uint64_t **[decompRNS];
        qInvModqpWithQP[lvlP] = new uint64_t ***[decompRNS];
        QModqpWithQP[lvlP] = new uint64_t **[decompRNS];
        for (int i = 0; i < decompRNS; i++) {
            qInvModqWithQP[lvlP][i] = new uint64_t *[xnbPi[i] - 1];
            qInvModqpWithQP[lvlP][i] = new uint64_t **[xnbPi[i] - 1];
            QModqpWithQP[lvlP][i] = new uint64_t *[xnbPi[i] - 1];
            for (int j = 0; j < xnbPi[i] - 1; j++) {
                QiWithQP.resize(j + 2);
                PiWithQP.resize(L + nbPi);
                for (int k = 0; k < j + 2; k++) {
                    QiWithQP[k] = qVec[i * nbPi + k];
                }
                for (int k = 0; k < L; k++) {
                    PiWithQP[k] = qVec[k];
                }
                for (int k = L; k < L + nbPi; k++) {
                    PiWithQP[k] = pVec[k - L];
                }
                qInvModqWithQP[lvlP][i][j] = new uint64_t[QiWithQP.size()];
                qInvModqpWithQP[lvlP][i][j] = new uint64_t *[PiWithQP.size()];
                for (int m = 0; m < PiWithQP.size(); m++) {
                    qInvModqpWithQP[lvlP][i][j][m] = new uint64_t[QiWithQP.size()];
                }
                for (int m = 0; m < QiWithQP.size(); m++) {
                    uint64_t qiStar = 1;
                    for (int n = 0; n < QiWithQP.size(); n++) {
                        if (n != m) {
                            uint64_t temp = QiWithQP[n] % QiWithQP[m];
                            MulMod(qiStar, qiStar, temp, QiWithQP[m]);
                        }
                    }
                    qInvModqWithQP[lvlP][i][j][m] = InvMod(qiStar, QiWithQP[m]);
                    for (int n = 0; n < PiWithQP.size(); n++) {
                        uint64_t qiStar = 1;
                        for (int u = 0; u < QiWithQP.size(); u++) {
                            if (u != m) {
                                uint64_t temp = QiWithQP[u] % PiWithQP[n];
                                MulMod(qiStar, qiStar, temp, PiWithQP[n]);
                            }
                        }
                        qInvModqpWithQP[lvlP][i][j][n][m] = qiStar;
                    }
                }
                QModqpWithQP[lvlP][i][j] = new uint64_t[PiWithQP.size()];
                for (int n = 0; n < PiWithQP.size(); n++) {
                    uint64_t QmodPi = 1;
                    for (int m = 0; m < QiWithQP.size(); m++) {
                        MulMod(QmodPi, QmodPi, QiWithQP[m], PiWithQP[n]);
                    }
                    QModqpWithQP[lvlP][i][j][n] = QmodPi;
                }
            }
        }
    }
    pSliceHatModp = new uint64_t *[K];
    for (long k = 0; k < K; k++) {
        pSliceHatModp[k] = new uint64_t[k + 1];
        for (long j = 0; j < k + 1; j++) {
            pSliceHatModp[k][j] = 1;
            for (long i = 0; i < j; i++) {
                uint64_t temp = pVec[i] % pVec[j];
                MulMod(pSliceHatModp[k][j], pSliceHatModp[k][j], temp, pVec[j]);
            }
            for (long i = j + 1; i < k + 1; i++) {
                uint64_t temp = pVec[i] % pVec[j];
                MulMod(pSliceHatModp[k][j], pSliceHatModp[k][j], temp, pVec[j]);
            }
        }
    }
    pSliceHatInvModp = new uint64_t *[K];
    for (long k = 0; k < K; k++) {
        pSliceHatInvModp[k] = new uint64_t[k + 1];
        for (long j = 0; j < k + 1; j++) {
            pSliceHatInvModp[k][j] = InvMod(pSliceHatModp[k][j], pVec[j]);
        }
    }
    pSliceHatModq = new uint64_t **[K];
    for (long k = 0; k < K; k++) {
        pSliceHatModq[k] = new uint64_t *[k + 1];
        for (long j = 0; j < k + 1; j++) {
            pSliceHatModq[k][j] = new uint64_t[L];
            for (long i = 0; i < L; i++) {
                pSliceHatModq[k][j][i] = 1;
                for (long s = 0; s < j; s++) {
                    uint64_t temp = pVec[s] % qVec[i];
                    MulMod(pSliceHatModq[k][j][i], pSliceHatModq[k][j][i], temp, qVec[i]);
                }
                for (long s = j + 1; s < k + 1; s++) {
                    uint64_t temp = pVec[s] % qVec[i];
                    MulMod(pSliceHatModq[k][j][i], pSliceHatModq[k][j][i], temp, qVec[i]);
                }
            }
        }
    }
    PSliceModq = new uint64_t *[K];
    for (long j = 0; j < K; j++) {
        PSliceModq[j] = new uint64_t[L];
        for (long i = 0; i < L; i++) {
            PSliceModq[j][i] = 1;
            for (long k = 0; k < j + 1; k++) {
                uint64_t temp = pVec[k] % qVec[i];
                MulMod(PSliceModq[j][i], PSliceModq[j][i], temp, qVec[i]);
            }
        }
    }
    PSliceInvModq = new uint64_t *[K];
    for (long j = 0; j < K; j++) {
        PSliceInvModq[j] = new uint64_t[L];
        for (long i = 0; i < L; i++) {
            PSliceInvModq[j][i] = InvMod(PSliceModq[j][i], qVec[i]);
        }
    }
    rescaleConstants = new uint64_t *[L - 1];
    for (int j = L - 1; j > 0; j--) {
        uint64_t qj = qVec[j];
        rescaleConstants[j - 1] = new uint64_t[j];
        for (int i = 0; i < j; i++) {
            uint64_t qi = qVec[i];
            rescaleConstants[j - 1][i] = qi - InvMod(qj, qi);
        }
    }
}

Context Context::InitContext(long logN, int lambda, int dNum)
{
    long L = 15;
    long K = 1;
    long T = 1;
    return Context(logN, L, K, T, lambda, dNum);
}

/*
 * @brief :对多项式在模数Q下进行数论变换 (NTT),并将结果存储到输出数组中。
 *
 * @param res :输出数组
 * @param a :输入的多项式数组
 * @param index :模数层级的索引
 */
void Context::QiNTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    QiNTTAndEqual(res, index);
}

/*
 * @brief :对多项式在模数P下进行数论变换 (NTT),并将结果存储到输出数组中。
 *
 * @param res :输出数组
 * @param a :输入的多项式数组
 * @param index :模数层级的索引
 */
void Context::PiNTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    PiNTTAndEqual(res, index);
}

/*
 * @brief :对多项式在模数PQ下进行数论变换 (NTT),并将结果存储到输出数组中。
 *
 * @param res :输出数组
 * @param a :输入的多项式数组
 * @param l: 模数Q层级数量。
 * @param k: 模数P层级数量。
 */
void Context::NTT(uint64_t *res, uint64_t *a, long l, long k)
{
#pragma omp parallel for 
    for (long index = 0; index < l; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        QiNTT(resi, ai, index);
    }
    for (long index = l; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        PiNTT(resi, ai, index - l);
    }
}

/*
 * @brief :在模数Q指定模数下,对多项式数组执行数论变换
 * (NTT),并将结果存储回原数组中
 *
 * @param a :输入/输出数组,包含待变换的多项式数据,结果会直接更新到该数组中。
 * @param index :模数索引
 */
void Context::QiNTTAndEqual(uint64_t *a, long index)
{
    long t = N;
    long logt1 = logN;
    uint64_t q = qVec[index];
    int gap = N / 8;
    for (long m = 1; m < gap; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = qRootPows[index][m + i];
            uint64_t wShoup = qRootPows_shoup[index][m + i];
            uint64x2_t q2 = vdupq_n_u64(q);
            for (long j = j1; j < j2; j += 4) {
                unsigned __int128 vMul[4];
                uint64_t hi[4];
                uint64_t V[4];
                for (int k = 0; k < 4; k++) {
                    vMul[k] = static_cast<unsigned __int128>(a[j + t + k]) * wShoup;
                    hi[k] = vMul[k] >> 64;
                    V[k] = a[j + t + k] * W - hi[k] * q;
                    if (V[k] >= q) {
                        V[k] -= q;
                    }
                }
                uint64x2_t vVec1 = vcombine_u64(vdup_n_u64(V[0]), vdup_n_u64(V[1]));
                uint64x2_t vVec2 = vcombine_u64(vdup_n_u64(V[2]), vdup_n_u64(V[3]));
                for (int k = 0; k < 2; k++) {
                    uint64x2_t ajtVec = vld1q_u64(&a[j + t + 2 * k]);
                    uint64x2_t ajVec = vld1q_u64(&a[j + 2 * k]);
                    uint64x2_t ajPlusQMinusV = vsubq_u64(vaddq_u64(ajVec, q2), k == 0 ? vVec1 : vVec2);
                    uint64x2_t lessMask = vcltq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t ajtNew = vbslq_u64(lessMask, ajPlusQMinusV, vsubq_u64(ajVec, k == 0 ? vVec1 : vVec2));
                    vst1q_u64(&a[j + t + 2 * k], ajtNew);
                    ajVec = vaddq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t gtMask = vcgtq_u64(ajVec, q2);
                    ajVec = vbslq_u64(gtMask, vsubq_u64(ajVec, q2), ajVec);
                    vst1q_u64(&a[j + 2 * k], ajVec);
                }
            }
        }
        logt1 -= 1;
    }
    for (long m = gap; m < N; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = qRootPows[index][m + i];
            uint64_t wShoup = qRootPows_shoup[index][m + i];
            for (long j = j1; j < j2; j++) {
                unsigned __int128 v_w_mul = static_cast<unsigned __int128>(a[j + t]) * wShoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = a[j + t] * W - hi * q;
                if (V >= q) {
                    V -= q;
                }
                a[j + t] = a[j] < V ? a[j] + q - V : a[j] - V;
                a[j] += V;
                if (a[j] >= q) {
                    a[j] -= q;
                }
            }
        }
        logt1 -= 1;
    }
}


/*
 * @brief :在模数P指定模数下,对多项式数组执行数论变换
 * (NTT),并将结果存储回原数组中
 *
 * @param a :输入/输出数组,包含待变换的多项式数据,结果会直接更新到该数组中。
 * @param index :模数索引
 */
void Context::PiNTTAndEqual(uint64_t *a, long index)
{
    long t = N;
    long logt1 = logN;
    uint64_t p = pVec[index];
    int gap = N / 8;
    for (long m = 1; m < gap; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = pRootPows[index][m + i];
            uint64_t wShoup = pRootPows_shoup[index][m + i];
            uint64x2_t q2 = vdupq_n_u64(p);
            for (long j = j1; j < j2; j += 4) {
                unsigned __int128 vMul[4];
                uint64_t hi[4];
                uint64_t V[4];
                for (int k = 0; k < 4; k++) {
                    vMul[k] = static_cast<unsigned __int128>(a[j + t + k]) * wShoup;
                    hi[k] = vMul[k] >> 64;
                    V[k] = a[j + t + k] * W - hi[k] * p;
                    if (V[k] >= p) {
                        V[k] -= p;
                    }
                }
                uint64x2_t vVec1 = vcombine_u64(vdup_n_u64(V[0]), vdup_n_u64(V[1]));
                uint64x2_t vVec2 = vcombine_u64(vdup_n_u64(V[2]), vdup_n_u64(V[3]));
                for (int k = 0; k < 2; k++) {
                    uint64x2_t ajtVec = vld1q_u64(&a[j + t + 2 * k]);
                    uint64x2_t ajVec = vld1q_u64(&a[j + 2 * k]);
                    uint64x2_t ajPlusQMinusV = vsubq_u64(vaddq_u64(ajVec, q2), k == 0 ? vVec1 : vVec2);
                    uint64x2_t lessMask = vcltq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t ajtNew = vbslq_u64(lessMask, ajPlusQMinusV, vsubq_u64(ajVec, k == 0 ? vVec1 : vVec2));
                    vst1q_u64(&a[j + t + 2 * k], ajtNew);
                    ajVec = vaddq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t gtMask = vcgtq_u64(ajVec, q2);
                    ajVec = vbslq_u64(gtMask, vsubq_u64(ajVec, q2), ajVec);
                    vst1q_u64(&a[j + 2 * k], ajVec);
                }
            }
        }
        logt1 -= 1;
    }
    for (long m = gap; m < N; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = pRootPows[index][m + i];
            uint64_t wShoup = pRootPows_shoup[index][m + i];
            for (long j = j1; j < j2; j++) {
                unsigned __int128 v_w_mul = static_cast<unsigned __int128>(a[j + t]) * wShoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = a[j + t] * W - hi * p;
                if (V >= p) {
                    V -= p;
                }
                a[j + t] = a[j] < V ? a[j] + p - V : a[j] - V;
                a[j] += V;
                if (a[j] >= p) {
                    a[j] -= p;
                }
            }
        }
        logt1 -= 1;
    }
}

/*
 * @brief :在模数T指定模数下,对多项式数组执行数论变换
 * (NTT),并将结果存储回原数组中
 *
 * @param a :输入/输出数组,包含待变换的多项式数据,结果会直接更新到该数组中。
 * @param index :模数索引
 */
void Context::TiNTTAndEqual(uint64_t *a, long index)
{
    long t = N;
    long logt1 = logN;
    uint64_t ti = tVec[index];
    int gap = N / 8;
    for (long m = 1; m < gap; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = tRootPows[index][m + i];
            uint64_t wShoup = tRootPows_shoup[index][m + i];
            uint64x2_t q2 = vdupq_n_u64(ti);
            for (long j = j1; j < j2; j += 4) {
                unsigned __int128 vMul[4];
                uint64_t hi[4];
                uint64_t V[4];
                for (int k = 0; k < 4; k++) {
                    vMul[k] = static_cast<unsigned __int128>(a[j + t + k]) * wShoup;
                    hi[k] = vMul[k] >> 64;
                    V[k] = a[j + t + k] * W - hi[k] * ti;
                    if (V[k] >= ti) {
                        V[k] -= ti;
                    }
                }
                uint64x2_t vVec1 = vcombine_u64(vdup_n_u64(V[0]), vdup_n_u64(V[1]));
                uint64x2_t vVec2 = vcombine_u64(vdup_n_u64(V[2]), vdup_n_u64(V[3]));
                for (int k = 0; k < 2; k++) {
                    uint64x2_t ajtVec = vld1q_u64(&a[j + t + 2 * k]);
                    uint64x2_t ajVec = vld1q_u64(&a[j + 2 * k]);
                    uint64x2_t ajPlusQMinusV = vsubq_u64(vaddq_u64(ajVec, q2), k == 0 ? vVec1 : vVec2);
                    uint64x2_t lessMask = vcltq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t ajtNew = vbslq_u64(lessMask, ajPlusQMinusV, vsubq_u64(ajVec, k == 0 ? vVec1 : vVec2));
                    vst1q_u64(&a[j + t + 2 * k], ajtNew);
                    ajVec = vaddq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t gtMask = vcgtq_u64(ajVec, q2);
                    ajVec = vbslq_u64(gtMask, vsubq_u64(ajVec, q2), ajVec);
                    vst1q_u64(&a[j + 2 * k], ajVec);
                }
            }
        }
        logt1 -= 1;
    }
    for (long m = gap; m < N; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = tRootPows[index][m + i];
            uint64_t wShoup = tRootPows_shoup[index][m + i];
            for (long j = j1; j < j2; j++) {
                unsigned __int128 v_w_mul = static_cast<unsigned __int128>(a[j + t]) * wShoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = a[j + t] * W - hi * ti;
                if (V >= ti) {
                    V -= ti;
                }
                a[j + t] = a[j] < V ? a[j] + ti - V : a[j] - V;
                a[j] += V;
                if (a[j] >= ti) {
                    a[j] -= ti;
                }
            }
        }
        logt1 -= 1;
    }
}

/*
 * @brief :对输入多项式数组在模数PQ下执行原地数论变换 (NTT),并更新数组数据。
 *
 * @param a :输入/输出数组
 * @param l: 模数Q层级数量。
 * @param k: 模数P层级数量。
 */
void Context::NTTAndEqual(uint64_t *a, long l, long k)
{
#pragma omp parallel for
    for (long index = 0; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        if (index < l) {
            QiNTTAndEqual(ai, index);
        } else {
            PiNTTAndEqual(ai, index - l);
        }
    }
}


/*
 * @brief :对多项式在模数Q下进行逆数论变换 (INTT),并将结果存储到输出数组中。
 *
 * @param res :输出数组
 * @param a :输入的多项式数组
 * @param index :模数层级的索引
 */
void Context::QiINTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    QiINTTAndEqual(res, index);
}

/*
 * @brief :对多项式在模数P下进行逆数论变换 (INTT),并将结果存储到输出数组中。
 *
 * @param res :输出数组
 * @param a :输入的多项式数组
 * @param index :模数层级的索引
 */
void Context::PiINTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    PiINTTAndEqual(res, index);
}

/*
 * @brief :在模数T指定模数下,对多项式数组执行逆数论变换
 * (INTT),并将结果存储回原数组中
 *
 * @param a :输入/输出数组,包含待变换的多项式数据,结果会直接更新到该数组中。
 * @param index :模数索引
 */
void Context::TiINTTAndEqual(uint64_t *a, long index)
{
    uint64_t ti = tVec[index];
    uint64_t td = tdVec[index];
    uint64_t tInv = tInvVec[index];
    long t = 1;
    for (long m = N; m > 1; m >>= 1) {
        long j1 = 0;
        long hm = m >> 1;
        for (long i = 0; i < hm; i++) {
            long j2 = j1 + t - 1;
            uint64_t W = tRootPowsInv[index][hm + i];
            uint64_t W_shoup = tRootPowsInv_shoup[index][hm + i];
            for (long j = j1; j <= j2; j++) {
                uint64_t Ti = a[j] + td;
                Ti -= a[j + t];
                a[j] += a[j + t];
                if (a[j] >= td)
                    a[j] -= td;
                unsigned __int128 v_w_mul = static_cast<unsigned __int128>(Ti) * W_shoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = Ti * W - hi * ti;
                if (V >= ti)
                    V -= ti;
                a[j + t] = V;
            }
            j1 += (t << 1);
        }
        t <<= 1;
    }
    uint64_t NScale = NInvModt[index];
    uint64_t NScale_shoup = NInvModt_shoup[index];
    uint64x2_t ti_vec = vdupq_n_u64(ti);
    uint64x2_t NScale_vec = vdupq_n_u64(NScale);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(&a[i]);
        uint64x2_t Ti_vec = vsubq_u64(a_vec, ti_vec);
        uint64x2_t mask = vcgeq_u64(a_vec, ti_vec);
        Ti_vec = vbslq_u64(mask, Ti_vec, a_vec);
        uint64_t Ti1 = vgetq_lane_u64(Ti_vec, 0);
        unsigned __int128 v_w_mul1 = static_cast<unsigned __int128>(Ti1) * NScale_shoup;
        uint64_t hi1 = v_w_mul1 >> 64;
        uint64_t V1 = Ti1 * NScale - hi1 * ti;
        if (V1 >= ti) {
            V1 -= ti;
        }
        uint64_t Ti2 = vgetq_lane_u64(Ti_vec, 1);
        unsigned __int128 v_w_mul2 = static_cast<unsigned __int128>(Ti2) * NScale_shoup;
        uint64_t hi2 = v_w_mul2 >> 64;
        uint64_t V2 = Ti2 * NScale - hi2 * ti;
        if (V2 >= ti) {
            V2 -= ti;
        }
        uint64x2_t V_vec = { V1, V2 };
        vst1q_u64(&a[i], V_vec);
    }
}

/*
 * @brief :在模数Q指定模数下,对多项式数组执行逆数论变换
 * (INTT),并将结果存储回原数组中
 *
 * @param a :输入/输出数组,包含待变换的多项式数据,结果会直接更新到该数组中。
 * @param index :模数索引
 */
void Context::QiINTTAndEqual(uint64_t *a, long index)
{
    uint64_t q = qVec[index];
    uint64_t qd = qdVec[index];
    uint64_t qInv = qInvVec[index];
    long t = 1;
    for (long m = N; m > 1; m >>= 1) {
        long j1 = 0;
        long hm = m >> 1;
        for (long i = 0; i < hm; i++) {
            long j2 = j1 + t - 1;
            uint64_t W = qRootPowsInv[index][hm + i];
            uint64_t W_shoup = qRootPowsInv_shoup[index][hm + i];
            for (long j = j1; j <= j2; j++) {
                uint64_t T = a[j] + qd;
                T -= a[j + t];
                a[j] += a[j + t];
                if (a[j] >= qd)
                    a[j] -= qd;
                unsigned __int128 v_w_mul = static_cast<unsigned __int128>(T) * W_shoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = T * W - hi * q;
                if (V >= q)
                    V -= q;
                a[j + t] = V;
            }
            j1 += (t << 1);
        }
        t <<= 1;
    }
    uint64_t NScale = NInvModq[index];
    uint64_t NScale_shoup = NInvModq_shoup[index];
    uint64x2_t qi_vec = vdupq_n_u64(q);
    uint64x2_t NScale_vec = vdupq_n_u64(NScale);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(&a[i]);
        uint64x2_t Ti_vec = vsubq_u64(a_vec, qi_vec);
        uint64x2_t mask = vcgeq_u64(a_vec, qi_vec);
        Ti_vec = vbslq_u64(mask, Ti_vec, a_vec);
        uint64_t Ti1 = vgetq_lane_u64(Ti_vec, 0);
        unsigned __int128 v_w_mul1 = static_cast<unsigned __int128>(Ti1) * NScale_shoup;
        uint64_t hi1 = v_w_mul1 >> 64;
        uint64_t V1 = Ti1 * NScale - hi1 * q;
        if (V1 >= q) {
            V1 -= q;
        }
        uint64_t Ti2 = vgetq_lane_u64(Ti_vec, 1);
        unsigned __int128 v_w_mul2 = static_cast<unsigned __int128>(Ti2) * NScale_shoup;
        uint64_t hi2 = v_w_mul2 >> 64;
        uint64_t V2 = Ti2 * NScale - hi2 * q;
        if (V2 >= q) {
            V2 -= q;
        }
        uint64x2_t V_vec = { V1, V2 };
        vst1q_u64(&a[i], V_vec);
    }
}

/*
 * @brief :在模数P指定模数下,对多项式数组执行逆数论变换
 * (INTT),并将结果存储回原数组中
 *
 * @param a :输入/输出数组,包含待变换的多项式数据,结果会直接更新到该数组中。
 * @param index :模数索引
 */
void Context::PiINTTAndEqual(uint64_t *a, long index)
{
    uint64_t p = pVec[index];
    uint64_t pd = pdVec[index];
    uint64_t pInv = pInvVec[index];
    long t = 1;
    for (long m = N; m > 1; m >>= 1) {
        long j1 = 0;
        long hm = m >> 1;
        for (long i = 0; i < hm; i++) {
            long j2 = j1 + t - 1;
            uint64_t W = pRootPowsInv[index][hm + i];
            uint64_t W_shoup = pRootPowsInv_shoup[index][hm + i];
            for (long j = j1; j <= j2; j++) {
                uint64_t T = a[j] + pd;
                T -= a[j + t];
                a[j] += a[j + t];
                if (a[j] >= pd)
                    a[j] -= pd;
                unsigned __int128 v_w_mul = static_cast<unsigned __int128>(T) * W_shoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = T * W - hi * p;
                if (V >= p)
                    V -= p;
                a[j + t] = V;
            }
            j1 += (t << 1);
        }
        t <<= 1;
    }
    uint64_t NScale = NInvModp[index];
    uint64_t NScale_shoup = NInvModp_shoup[index];
    uint64x2_t pi_vec = vdupq_n_u64(p);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(&a[i]);
        uint64x2_t Ti_vec = vsubq_u64(a_vec, pi_vec);
        uint64x2_t mask = vcgeq_u64(a_vec, pi_vec);
        Ti_vec = vbslq_u64(mask, Ti_vec, a_vec);
        uint64_t Ti1 = vgetq_lane_u64(Ti_vec, 0);
        unsigned __int128 v_w_mul1 = static_cast<unsigned __int128>(Ti1) * NScale_shoup;
        uint64_t hi1 = v_w_mul1 >> 64;
        uint64_t V1 = Ti1 * NScale - hi1 * p;
        if (V1 >= p) {
            V1 -= p;
        }
        uint64_t Ti2 = vgetq_lane_u64(Ti_vec, 1);
        unsigned __int128 v_w_mul2 = static_cast<unsigned __int128>(Ti2) * NScale_shoup;
        uint64_t hi2 = v_w_mul2 >> 64;
        uint64_t V2 = Ti2 * NScale - hi2 * p;
        if (V2 >= p) {
            V2 -= p;
        }
        uint64x2_t V_vec = { V1, V2 };
        vst1q_u64(&a[i], V_vec);
    }
}

/*
 * @brief : 对输入多项式数组执行原地逆数论变换 (INTT),并更新数组数据。
 *
 * @param a :输入/输出数组
 * @param l: 模数Q层级数量。
 * @param k: 模数P层级数量。
 */
void Context::INTTAndEqual(uint64_t *a, long l, long k)
{
#pragma omp parallel for
    for (long index = 0; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        if (index < l) {
            QiINTTAndEqual(ai, index);
        } else {
            PiINTTAndEqual(ai, index - l);
        }   
    }
}

/*
 * @brief :对多项式在模数PQ下进行逆数论变换 (NTT),并将结果存储到输出数组中。
 *
 * @param res :输出数组
 * @param a :输入的多项式数组
 * @param l: 模数Q层级数量。
 * @param k: 模数P层级数量。
 */
void Context::INTT(uint64_t *res, uint64_t *a, long l, long k)
{
#pragma omp parallel for
    for (long index = 0; index < l; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        QiINTT(resi, ai, index);
    }
    for (long index = l; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        PiINTT(resi, ai, index - l);
    }
}

/*
 * @brief 在模数Q指定模数下,对输入数组中的多项式系数进行取负操作。
 *
 * @param res: 输出数组,用于存储取负后的结果。
 * @param a: 输入数组,包含待取负的多项式系数。
 * @param index: 模数索引,指定当前取负操作所使用的模数。
 */
void Context::QiNegate(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        res[i] = (a[i] == 0) ? 0 : qVec[index] - a[i];
    }
}

/*
 * @brief 在模数P指定模数下,对输入数组中的多项式系数进行取负操作。
 *
 * @param res: 输出数组,用于存储取负后的结果。
 * @param a: 输入数组,包含待取负的多项式系数。
 * @param index: 模数索引,指定当前取负操作所使用的模数。
 */
void Context::PiNegate(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        res[i] = (a[i] == 0) ? 0 : pVec[index] - a[i];
    }
}

/*
 * @brief 在指定模数下,对输入多项式数组执行取负操作,并将结果存储到输出数组中。
 *
 * @param res: 输出数组,用于存储取负后的多项式系数。
 * @param a: 输入数组,包含待取负的多项式系数。
 * @param l: 处理的模数层级数量,前 `l` 层调用 `QiNegate` 进行取负操作。
 * @param k: 额外层级数量,后 `k` 层调用 `PiNegate` 进行取负操作。
 */
void Context::Negate(uint64_t *res, uint64_t *a, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiNegate(resi, ai, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiNegate(resi, ai, i - l);
    }
}

/*
 * @brief 在模数Q指定模数下,对输入数组中的多项式系数进行取负操作。
 *
 * @param a: 输入/输出数组。
 * @param index: 模数索引,指定当前取负操作所使用的模数。
 */
void Context::QiNegateAndEqual(uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        a[i] = (a[i] == 0) ? 0 : qVec[index] - a[i];
    }
}

/*
 * @brief 在模数P指定模数下,对输入数组中的多项式系数进行取负操作。
 *
 * @param a: 输入/输出数组。
 * @param index: 模数索引,指定当前取负操作所使用的模数。
 */
void Context::PiNegateAndEqual(uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        a[i] = (a[i] == 0) ? 0 : pVec[index] - a[i];
    }
}

/*
 * @brief 在指定模数下,对输入多项式数组执行取负操作,并将结果存储到输出数组中。
 *
 * @param a: 输入/输出数组。
 * @param l: 处理的模数层级数量,前 `l` 层调用 `QiNegate` 进行取负操作。
 * @param k: 额外层级数量,后 `k` 层调用 `PiNegate` 进行取负操作。
 */
void Context::NegateAndEqual(uint64_t *a, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        QiNegateAndEqual(ai, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        PiNegateAndEqual(ai, i - l);
    }
}

/*
 * @brief 在模数Q指定模数下,对输入多项式数组的每个元素加上一个常数,并存储结果。
 *
 * @param res: 输出数组,用于存储加上常数后的多项式系数。
 * @param a: 输入数组,包含待加常数的多项式系数。
 * @param c: 要加上的常数值。
 * @param index: 模数索引,用于从 `qVec` 中获取当前模数。
 */
void Context::QiAddConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        AddMod(res[i], a[i], c, qVec[index]);
    }
}

/*
 * @brief 在模数P指定模数下,对输入多项式数组的每个元素加上一个常数,并存储结果。
 *
 * @param res: 输出数组,用于存储加上常数后的多项式系数。
 * @param a: 输入数组,包含待加常数的多项式系数。
 * @param c: 要加上的常数值。
 * @param index: 模数索引,用于从 `qVec` 中获取当前模数。
 */
void Context::PiAddConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        AddMod(res[i], a[i], c, pVec[index]);
    }
}

/*
 * @brief 在指定模数域下,对输入多项式数组的每个元素加上一个常数,并存储结果。
 *
 * @param res: 输出数组,用于存储加常数后的多项式系数。
 * @param a: 输入数组,包含待加常数的多项式系数。
 * @param c: 要加上的常数值。
 * @param l: 处理的模数层级数量,前 `l` 层调用 `QiAddConst` 进行加法操作。
 * @param k: 额外层级数量,后 `k` 层调用 `PiAddConst` 进行加法操作。
 */
void Context::AddConst(uint64_t *res, uint64_t *a, uint64_t c, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiAddConst(resi, ai, c, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiAddConst(resi, ai, c, i - l);
    }
}

/*
 * @brief
 * 在指定模数下,对输入多项式数组的每个元素加上一个常数,结果直接更新到原数组中。
 *
 * @param a: 输入/输出数组,包含待加常数的多项式系数,结果将直接更新到此数组。
 * @param c: 要加上的常数值。
 * @param index: 模数索引,用于从 `qVec` 中获取当前的模数。
 */
void Context::QiAddConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        AddMod(a[i], a[i], c, qVec[index]);
    }
}

/*
 * @brief
 * 在指定额外模数下,对输入多项式数组的每个元素加上一个常数,结果直接更新到原数组中。
 *
 * @param a: 输入/输出数组,包含待加常数的多项式系数,结果将直接更新到此数组。
 * @param c: 要加上的常数值。
 * @param index: 额外模数索引,用于从 `pVec` 中获取当前的模数。
 */
void Context::PiAddConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    // #pragma omp parallel for
    for (long i = 0; i < N; i++) {
        AddMod(a[i], a[i], c, pVec[index]);
    }
}

/*
 * @brief 在指定模数域下，对输入多项式数组的每个元素加上一个常数，并存储结果。
 *
 * @param a： 输入/输出数组，包含待加常数的多项式系数，结果将直接更新到此数组。
 * @param c： 要加上的常数值。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiAddConstAndEqual`
 * 进行加法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiAddConstAndEqual` 进行加法操作。
 */
void Context::AddConstAndEqual(uint64_t *a, uint64_t c, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        QiAddConstAndEqual(ai, c, i);
    }

    for (long i = L; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        PiAddConstAndEqual(ai, c, i - l);
    }
}

/*
 * @brief
 * 在模数Q指定模数下，对输入多项式数组的每个元素减去一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储减去常数后的多项式系数。
 * @param a： 输入数组，包含待减常数的多项式系数。
 * @param c： 要减去的常数值。
 * @param index： 模数索引，用于从 `qVec` 中获取当前模数。
 */
void Context::QiSubConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], c, qVec[index]);
    }
}

/*
 * @brief
 * 在模数P指定模数下，对输入多项式数组的每个元素减去一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储减去常数后的多项式系数。
 * @param a： 输入数组，包含待减常数的多项式系数。
 * @param c： 要减去的常数值。
 * @param index： 模数索引，用于从 `pVec` 中获取当前模数。
 */
void Context::PiSubConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], c, pVec[index]);
    }
}

/*
 * @brief 在指定模数域下，对输入多项式数组的每个元素减去一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储减去常数后的多项式系数。
 * @param a： 输入数组，包含待减常数的多项式系数。
 * @param c： 要减去的常数值。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiSubConst` 进行减法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiSubConst` 进行减法操作。
 */
void Context::SubConst(uint64_t *res, uint64_t *a, uint64_t c, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiSubConst(resi, ai, c, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiSubConst(resi, ai, c, i - l);
    }
}

/*
 * @brief
 * 在指定模数下，对输入多项式数组的每个元素减去一个常数，结果直接更新到原数组中。
 *
 * @param a： 输入/输出数组，包含待减常数的多项式系数，结果将直接更新到此数组。
 * @param c： 要减去的常数值。
 * @param index： 模数索引，用于从 `qVec` 中获取当前的模数。
 */
void Context::QiSubConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(a[i], a[i], c, qVec[index]);
    }
}

/*
 * @brief
 * 在指定模数下，对输入多项式数组的每个元素减去一个常数，结果直接更新到原数组中。
 *
 * @param a： 输入/输出数组，包含待减常数的多项式系数，结果将直接更新到此数组。
 * @param c： 要减去的常数值。
 * @param index： 模数索引，用于从 `pVec` 中获取当前的模数。
 */
void Context::PiSubConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(a[i], a[i], c, pVec[index]);
    }
}

/*
 * @brief 在指定模数域下，对输入多项式数组的每个元素减去一个常数，并存储结果。
 *
 * @param a： 输入/输出数组，包含待减常数的多项式系数，结果将直接更新到此数组。
 * @param c： 要减去的常数值。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiSubConstAndEqual`
 * 进行减法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiSubConstAndEqual` 进行减法操作。
 */
void Context::SubConstAndEqual(uint64_t *a, uint64_t c, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        QiSubConstAndEqual(ai, c, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        PiSubConstAndEqual(ai, c, i - l);
    }
}

/*
 * @brief ：对多项式在模数Q下进行加法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiAdd(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t q_vec = vdupq_n_u64(qVec[index]);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, q_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, q_vec));
        vst1q_u64(res + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数P下进行加法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiAdd(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t p_vec = vdupq_n_u64(pVec[index]);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, p_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, p_vec));
        vst1q_u64(res + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数T下进行加法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `tVec` 中获取当前的模数
 */
void Context::TiAdd(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t t_vec = vdupq_n_u64(tVec[index]);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, t_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, t_vec));
        vst1q_u64(res + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数QP下进行加法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param b ：输入的多项式数组
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiAdd` 或 `QiAdd_NEON`
 * 进行加法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiAdd` 或 `PiAdd_NEON` 进行加法操作。
 */
void Context::Add(uint64_t *res, uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiAdd(resi, ai, bi, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiAdd(resi, ai, bi, i - l);
    }
}

/*
 * @brief ：对多项式在模数Q下进行加法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiAddAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t q_vec = vdupq_n_u64(qVec[index]);
#pragma omp parallel for
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, q_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, q_vec));
        vst1q_u64(a + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数P下进行加法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiAddAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t p_vec = vdupq_n_u64(pVec[index]);
#pragma omp parallel for
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, p_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, p_vec));
        vst1q_u64(a + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数T下进行加法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `tVec` 中获取当前的模数
 */
void Context::TiAddAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t t_vec = vdupq_n_u64(tVec[index]);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, t_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, t_vec));
        vst1q_u64(a + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数T下进行加法，并将结果存储到输出数组中。
 *
 * @param a ：输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param t： 处理的模数层级数量，调用 `TiAddAndEqual` 进行加法操作。
 */
void Context::TAddAndEqual(uint64_t *a, uint64_t *b, long t)
{
    for (long i = 0; i < t; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        TiAddAndEqual(ai, bi, i);
    }
}

/*
 * @brief ：对多项式在模数QP下进行加法，并将结果存储到输出数组中。
 *
 * @param a ：输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiAddAndEqual` 进行加法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiAddAndEqual` 进行加法操作。
 */
void Context::AddAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        if (i < l) {
            QiAddAndEqual(ai, bi, i);
        } else {
            PiAddAndEqual(ai, bi, i - l);
        }
    }
}

/*
 * @brief ：对多项式在模数Q下进行减法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiSub(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], b[i], qVec[index]);
    }
}

/*
 * @brief ：对多项式在模数P下进行减法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiSub(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], b[i], pVec[index]);
    }
}

/*
 * @brief ：对多项式在模数QP下进行减法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiSub` 进行减法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiSub` 进行减法操作。
 */
void Context::Sub(uint64_t *res, uint64_t *a, uint64_t *b, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiSub(resi, ai, bi, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiSub(resi, ai, bi, i - l);
    }
}

/*
 * @brief ：对多项式在模数PQ下进行减法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param k： 额外层级数量，前 `k` 层调用 `PiSub` 进行减法操作。
 * @param l： 处理的模数层级数量，后 `l` 层调用 `QiSub` 进行减法操作。
 */
void Context::TSub(uint64_t *res, uint64_t *a, uint64_t *b, long l, long k)
{
    for (long i = 0; i < k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiSub(resi, ai, bi, i);
    }
    for (long i = k; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiSub(resi, ai, bi, i - k);
    }
}

/*
 * @brief ：对多项式在模数Q下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiSubAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64_t r;
    uint64_t q = qVec[index];
    for (long i = 0; i < N; i++) {
        r = b[i] % q;
        a[i] = (a[i] + q - r) % q;
    }
}

/*
 * @brief ：对多项式在模数P下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiSubAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64_t r;
    uint64_t p = pVec[index];
    for (long i = 0; i < N; i++) {
        r = b[i] % p;
        a[i] = (a[i] + p - r) % p;
    }
}

/*
 * @brief ：对多项式在模数T下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param index ：模数层级的索引，用于从 `tVec` 中获取当前的模数
 */
void Context::TiSubAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64_t r;
    uint64_t t = tVec[index];
    for (long i = 0; i < N; i++) {
        r = b[i] % t;
        a[i] = (a[i] + t - r) % t;
    }
}

/*
 * @brief ：对多项式在模数QP下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiSubAndEqual` 进行减法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiSubAndEqual` 进行减法操作。
 */
void Context::SubAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        if (i < l) {
            QiSubAndEqual(ai, bi, i);
        } else {
            PiSubAndEqual(ai, bi, i - l);
        }
    }
}

/*
 * @brief ：对多项式在模数PQ下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。（被减数）
 * @param b ：输入的多项式数组（减数）
 * @param k： 额外层级数量，前 `k` 层调用 `PiSubAndEqual` 进行减法操作。
 * @param l： 处理的模数层级数量，后 `l` 层调用 `QiSubAndEqual` 进行减法操作。
 */
void Context::TSubAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
    for (long i = 0; i < k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        PiSubAndEqual(ai, bi, i);
    }
    for (long i = k; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        QiSubAndEqual(ai, bi, i - k);
    }
}

/*
 * @brief ：对多项式在模数Q下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入多项式数组（减数）
 * @param b ：输入/输出的多项式数组，结果将直接更新到此数组。（被减数）
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiSub2AndEqual(uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(b[i], a[i], b[i], qVec[index]);
    }
}

/*
 * @brief ：对多项式在模数P下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入多项式数组（减数）
 * @param b ：输入/输出的多项式数组，结果将直接更新到此数组。（被减数）
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiSub2AndEqual(uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(b[i], a[i], b[i], pVec[index]);
    }
}

/*
 * @brief ：对多项式在模数QP下进行减法，并将结果存储到输出数组中。
 *
 * @param a： 输入多项式数组（减数）
 * @param b ：输入/输出的多项式数组，结果将直接更新到此数组。（被减数）
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiSub2AndEqual` 进行减法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiSub2AndEqual` 进行减法操作。
 */
void Context::Sub2AndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        if (i < l) {
            QiSub2AndEqual(ai, bi, i);
        } else {
            PiSub2AndEqual(ai, bi, i - l);
        }
    }
}

/*
 * @brief
 * 在模数Q指定模数下，对输入多项式数组的每个元素乘上一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储乘上常数后的多项式系数。
 * @param a： 输入数组，包含待乘常数的多项式系数。
 * @param cnst： 要乘上的常数值。
 * @param index： 模数索引，用于从 `qVec` 中获取当前模数。
 */
void Context::QiMulConst(uint64_t *res, uint64_t *a, uint64_t cnst, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], cnst, qVec[index], twoqVec[index], qrVec_high[index], qrVec_low[index]);
    }
}

/*
 * @brief
 * 在模数P指定模数下，对输入多项式数组的每个元素乘上一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储乘上常数后的多项式系数。
 * @param a： 输入数组，包含待乘常数的多项式系数。
 * @param cnst： 要乘上的常数值。
 * @param index： 模数索引，用于从 `pVec` 中获取当前模数。
 */
void Context::PiMulConst(uint64_t *res, uint64_t *a, uint64_t cnst, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], cnst, pVec[index], twopVec[index], prVec_high[index], prVec_low[index]);
    }
}

/*
 * @brief 在指定模数域下，对输入多项式数组的每个元素乘上一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储乘常数后的多项式系数。
 * @param a： 输入数组，包含待乘常数的多项式系数。
 * @param cnst： 要乘上的常数值。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiMulConst` 进行乘法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiMulConst` 进行乘法操作。
 */
void Context::MulConst(uint64_t *res, uint64_t *a, uint64_t cnst, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiMulConst(resi, ai, cnst, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiMulConst(resi, ai, cnst, i - l);
    }
}

/*
 * @brief
 * 在指定模数Q下，对输入多项式数组的每个元素乘上一个常数，结果直接更新到原数组中。
 *
 * @param res：
 * 输入/输出数组，包含待乘常数的多项式系数，结果将直接更新到此数组。
 * @param cnst： 要乘上的常数值。
 * @param index： 模数索引，用于从 `qVec` 中获取当前的模数。
 */
void Context::QiMulConstAndEqual(uint64_t *res, uint64_t cnst, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], res[i], cnst, qVec[index], twoqVec[index], qrVec_high[index], qrVec_low[index]);
    }
}

/*
 * @brief
 * 在指定模数P下，对输入多项式数组的每个元素乘上一个常数，结果直接更新到原数组中。
 *
 * @param res：
 * 输入/输出数组，包含待乘常数的多项式系数，结果将直接更新到此数组。
 * @param cnst： 要乘上的常数值。
 * @param index： 模数索引，用于从 `pVec` 中获取当前的模数。
 */
void Context::PiMulConstAndEqual(uint64_t *res, uint64_t cnst, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], res[i], cnst, pVec[index], twopVec[index], prVec_high[index], prVec_low[index]);
    }
}

/*
 * @brief 在指定模数域下，对输入多项式数组的每个元素乘上一个常数，并存储结果。
 *
 * @param a： 输入/输出数组，包含待乘常数的多项式系数，结果将直接更新到此数组。
 * @param cnst： 要乘上的常数值。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiMulConstAndEqual`
 * 进行乘法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiMulConstAndEqual` 进行乘法操作。
 */
void Context::MulConstAndEqual(uint64_t *a, uint64_t cnst, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        if (i < l) {
            QiMulConstAndEqual(ai, cnst, i);
        } else {
            PiMulConstAndEqual(ai, cnst, i - l);
        }
    }
}

/*
 * @brief ：对多项式在模数Q下进行乘法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiMul(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    #pragma omp parallel for
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], b[i], qVec[index], twoqVec[index], qrVec_high[index], qrVec_low[index]);
    }
}

/*
 * @brief ：对多项式在模数P下进行乘法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiMul(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t p_vec = vdupq_n_u64(pVec[index]);
    uint64x2_t twop_vec = vdupq_n_u64(twopVec[index]);
    uint64x2_t pr_high_vec = vdupq_n_u64(prVec_high[index]);
    uint64x2_t pr_low_vec = vdupq_n_u64(prVec_low[index]);
#pragma omp parallel for
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t res_vec;
        MulModBarrett_NEON(&res_vec, a_vec, b_vec, p_vec, twop_vec, pr_high_vec, pr_low_vec);
        vst1q_u64(res + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数QP下进行乘法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param b ：输入的多项式数组
 * @param l： 处理的模数层级数量，前 `l` 层在模Q上进行乘法操作。
 * @param k： 额外层级数量，后 `k` 层在模P上进行乘法操作。
 */
void Context::Mul(uint64_t *res, uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiMul(resi, ai, bi, i);   
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiMul(resi, ai, bi, i - l);       
    }
}

/*
 * @brief ：对多项式在模数QP下进行乘法，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param key ：输入的多项式数组
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiMul` 进行乘法操作。
 */
void Context::MulKey(uint64_t *res, uint64_t *a, uint64_t *key, long l)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *keyi = key + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiMul(resi, ai, keyi, i);
    }
    for (long i = l; i < l + K; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *keyi = key + ((i - l + L) << logN);
        uint64_t *resi = res + (i << logN);
        PiMul(resi, ai, keyi, i - l);
    }
}

/*
 * @brief ：对多项式在模数Q下进行乘法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiMulAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t q_vec = vdupq_n_u64(qVec[index]);
    uint64x2_t twoq_vec = vdupq_n_u64(twoqVec[index]);
    uint64x2_t qr_high_vec = vdupq_n_u64(qrVec_high[index]);
    uint64x2_t qr_low_vec = vdupq_n_u64(qrVec_low[index]);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t res_vec;
        MulModBarrett_NEON(&res_vec, a_vec, b_vec, q_vec, twoq_vec, qr_high_vec, qr_low_vec);
        vst1q_u64(a + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数P下进行乘法，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiMulAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t p_vec = vdupq_n_u64(pVec[index]);
    uint64x2_t twop_vec = vdupq_n_u64(twopVec[index]);
    uint64x2_t pr_high_vec = vdupq_n_u64(prVec_high[index]);
    uint64x2_t pr_low_vec = vdupq_n_u64(prVec_low[index]);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t res_vec;
        MulModBarrett_NEON(&res_vec, a_vec, b_vec, p_vec, twop_vec, pr_high_vec, pr_low_vec);
        vst1q_u64(a + i, res_vec);
    }
}

/*
 * @brief ：对多项式在模数QP下进行乘法，并将结果存储到输出数组中。
 *
 * @param a ：输入/输出多项式数组，结果将直接更新到此数组。
 * @param b ：输入的多项式数组
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiMulAndEqual` 进行乘法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiMulAndEqual` 进行乘法操作。
 */
void Context::MulAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        if (i < l) {
            QiMulAndEqual(ai, bi, i);
        } else {
            PiMulAndEqual(ai, bi, i - l);
        }
    }
}
/*
 * @brief ：对多项式在模数Q下进行平方，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiSquare(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], a[i], qVec[index], twoqVec[index], qrVec_high[index], qrVec_low[index]);
    }
}

/*
 * @brief ：对多项式在模数P下进行平方，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiSquare(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], a[i], pVec[index], twopVec[index], prVec_high[index], prVec_low[index]);
    }
}

/*
 * @brief ：对多项式在模数QP下进行平方，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiSquare` 进行平方操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiSquare` 进行平方操作
 */
void Context::Square(uint64_t *res, uint64_t *a, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiSquare(resi, ai, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiSquare(resi, ai, i - l);
    }
}

/*
 * @brief ：对多项式在模数Q下进行平方，并将结果存储到输出数组中。
 *
 * @param a ：输入/输出的多项式数组，将更新结果直接写入到此数组中
 * @param index ：模数层级的索引，用于从 `qVec` 中获取当前的模数
 */
void Context::QiSquareAndEqual(uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(a[i], a[i], a[i], qVec[index], twoqVec[index], qrVec_high[index], qrVec_low[index]);
    }
}

/*
 * @brief ：对多项式在模数P下进行平方，并将结果存储到输出数组中。
 *
 * @param a ：输入/输出的多项式数组，将更新结果直接写入到此数组中
 * @param index ：模数层级的索引，用于从 `pVec` 中获取当前的模数
 */
void Context::PiSquareAndEqual(uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(a[i], a[i], a[i], pVec[index], twopVec[index], prVec_high[index], prVec_low[index]);
    }
}

/*
 * @brief ：对多项式在模数QP下进行平方，并将结果存储到输出数组中。
 *
 * @param a ：输入/输出的多项式数组，将更新结果直接写入到此数组中
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiSquare` 进行平方操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiSquare` 进行平方操作
 */
void Context::SquareAndEqual(uint64_t *a, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        QiSquare(ai, ai, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        PiSquare(ai, ai, i - l);
    }
}

/* *
 * @brief 对密文进行重缩放操作并更新密文。
 *
 * @param a 密文对象，将被更新。
 * @param l 重缩放的层数。
 */
void Context::ReScaleAndEqual(uint64_t *&a, long l)
{
    uint64_t *ra = new uint64_t[(l - 1) << logN];
    uint64_t *al = a + ((l - 1) << logN);
    QiINTTAndEqual(al, l - 1);
    uint64_t pj = qVec[l - 1];
    uint64_t pHalf = (pj - 1) >> 1;
    QiAddConstAndEqual(al, pHalf, l - 1);
#pragma omp parallel for
    for (long i = 0; i < l - 1; i++) {
        uint64_t qi = qVec[i];
        uint64_t *rai = ra + (i << logN);
        uint64_t *ai = a + (i << logN);
        uint64_t tmp = 0;
        ModBarrett(tmp, pHalf, qVec[i], twoqVec[i], qrVec_high[i], qrVec_low[i]);
        QiAddConst(rai, al, qi - tmp, i);
        QiNTTAndEqual(rai, i);
        for (long n = 0; n < N; n++) {
            rai[n] = (qVec[i] << 1) - ai[n] + rai[n];
            MulModBarrett(rai[n], rai[n], rescaleConstants[l - 2][i], qVec[i], twoqVec[i], qrVec_high[i], qrVec_low[i]);
        }
    }
    delete[] a;
    a = ra;
}

/*
 * @brief :生成服从高斯分布的多项式系数样本
 *
 * @param res: 输出数组,用于存储生成的高斯分布样本。
 * @param l: 处理的模数层级数量,前 `l` 层使用标准模数向量 `qVec`。
 * @param k: 额外层级数量,后 `k` 层使用额外模数向量 `pVec`。
 */
void Context::SampleGauss(uint64_t *res, long l, long k)
{
    // 使用 C++11 <random> 库的随机数生成器和分布
    // 创建静态随机数生成器，确保只初始化一次
    static std::mt19937_64 gen(std::random_device{}());

    // 定义均匀分布，范围 (0, 1)
    static std::uniform_real_distribution<double> dis(0.0, 1.0);
#pragma omp parallel for
    for (long i = 0; i < N; i += 2) {
        // 生成两个独立的均匀分布随机数，范围 (0, 1)
        double r1 = dis(gen);
        double r2 = dis(gen);

        // 使用 Box-Muller 变换将均匀分布转换为标准正态分布
        double theta = 2.0 * M_PI * r1;
        double rr = sqrt(-2.0 * log(r2)) * sigma;

        // 生成两个高斯随机数，并四舍五入到最近的整数
        long g1 = static_cast<long>(std::floor(rr * std::cos(theta) + 0.5));
        long g2 = static_cast<long>(std::floor(rr * std::sin(theta) + 0.5));

        // 将生成的高斯随机数应用到结果数组
        
        for (long j = 0; j < l; j++) {
            uint64_t *resj = res + (j << logN);
            resj[i] = g1 >= 0 ? static_cast<uint64_t>(g1) : qVec[j] + g1;
            resj[i + 1] = g2 >= 0 ? static_cast<uint64_t>(g2) : qVec[j] + g2;
        }
        for (long j = 0; j < k; j++) {
            uint64_t *resj = res + ((j + l) << logN);
            resj[i] = g1 >= 0 ? static_cast<uint64_t>(g1) : pVec[j] + g1;
            resj[i + 1] = g2 >= 0 ? static_cast<uint64_t>(g2) : pVec[j] + g2;
        }
    }
}

/*
 * @brief :生成服从 {0, 1, -1} 分布的多项式系数样本
 *
 * @param res: 输出数组,用于存储生成的 {0, +1, -1} 分布样本。
 * @param s: 未使用的参数（预留接口一致性,可能为以后扩展所用）。
 * @param l: 处理的模数层级数量,前 `l` 层使用标准模数向量 `qVec`。
 * @param k: 额外层级数量,后 `k` 层使用额外模数向量 `pVec`。
 */
void Context::SampleZO(uint64_t *res, long l, long k)
{
    // 使用 C++11 <random> 库的随机数生成器和分布
    // 创建静态随机数生成器，确保只初始化一次
    static std::mt19937_64 gen(std::random_device{}());

    // 定义分布
    // 对于 zo = {0, 1, -1}，每个值的概率均为1/3
    std::uniform_int_distribution<int> dist_zo(0, 2); // 0,1,2 分别映射到 0,1,-1

    for (long i = 0; i < N; i++) {
        int rand_val = dist_zo(gen);
        long zo;
        switch (rand_val) {
            case 0:
                zo = 0;
                break;
            case 1:
                zo = 1;
                break;
            case 2:
                zo = -1;
                break;
            default:
                zo = 0; // 理论上不会到达
        }

        if (std::abs(zo) > 1)
            std::cerr << "error in TUG !!!!!!!\n";

        for (long j = 0; j < l; j++) {
            uint64_t *resj = res + (j << logN);
            resj[i] = zo >= 0 ? static_cast<uint64_t>(zo) : qVec[j] + zo;
        }
        for (long j = 0; j < k; j++) {
            uint64_t *resj = res + ((j + l) << logN);
            resj[i] = zo >= 0 ? static_cast<uint64_t>(zo) : pVec[j] + zo;
        }
    }
}


/*
 * @brief :生成均匀分布的随机多项式系数样本。
 *
 * @param res: 输出数组,用于存储生成的均匀分布样本。
 * @param l: 处理的模数层级数量,前 `l` 层使用标准模数向量 `qVec`。
 * @param k: 额外层级数量,后 `k` 层使用额外模数向量 `pVec`。
 */
void Context::SampleUniform(uint64_t *res, long l, long k)
{
    // 使用 C++11 <random> 库的随机数生成器和分布
    // 创建静态随机数生成器，确保只初始化一次
    static std::mt19937_64 gen(std::random_device{}());

    // 定义均匀实数分布，范围 [0.0, 1.0)
    std::uniform_real_distribution<double> dis(0.0, 1.0);

    // 生成均匀分布的随机数并赋值
#pragma omp parallel for
    for (long j = 0; j < l; j++) {
        uint64_t *resj = res + (j << logN);
        for (long n = 0; n < N; n++) {
            double rand_real = dis(gen);
            resj[n] = static_cast<uint64_t>(std::floor(rand_real * qVec[j]));
        }
    }
    for (long j = 0; j < k; j++) {
        uint64_t *resj = res + ((j + l) << logN);
        for (long n = 0; n < N; n++) {
            double rand_real = dis(gen);
            resj[n] = static_cast<uint64_t>(std::floor(rand_real * pVec[j]));
        }
    }
}

/*
 * @brief :计算 (a^k) % mod 的快速幂结果,使用二进制快速幂算法。
 *
 * @param a: 底数,参与指数幂运算的数值。
 * @param k: 指数,表示幂次。
 * @param mod: 模数,用于取模运算,确保结果在 [0, mod) 范围内。
 *
 * @return uint64_t: (a^k) % mod 的计算结果。
 */
uint64_t Context::FastPow(uint64_t a, uint64_t k, uint64_t mod)
{
    uint64_t base = 1;
    __uint128_t b = a;
    if (mod == 0) {
        return 0;
    }
    while (k) {
        if (k & 1)
            base = (base * b) % mod;
        b = (b * b) % mod;
        k >>= 1;
    }
    return base % mod;
}

/*
 * @brief 把整数按二进制形式翻转
 *
 * @param index 输入要翻转的整数
 * @param bitLen 输入整数的二进制位数
 * @return uint64_t 返回翻转后的整数
 */
uint64_t Context::BitReverse64(uint64_t index, uint64_t bitLen)
{
    uint64_t reversed = 0;
    for (uint64_t i = 0; i < bitLen; i++) {
        reversed = (reversed << 1) | (index & 1);
        index >>= 1;
    }
    return reversed;
}

/*
 * @brief :生成 NTT 索引重排数组,基于输入的伽罗瓦元素进行重排。
 *
 * @param galEl :输入的伽罗瓦元素,用于计算索引重排的映射关系。
 * @param ntt_index :输出数组,用于存储重排后的 NTT 索引。
 */
void Context::PermuteNTTIndex(uint64_t galEl, uint64_t *ntt_index)
{
    uint64_t tmp1, tmp2;
    uint64_t logNthRoot = log2(M) - 1;
    uint64_t mask = M - 1;
    for (int i = 0; i < N; i++) {
        tmp1 = 2 * BitReverse64(i, logNthRoot) + 1;
        tmp2 = ((galEl * tmp1 & mask) - 1) >> 1;
        ntt_index[i] = BitReverse64(tmp2, logNthRoot);
    }
}

/*
 * @brief :使用给定的重排索引数组,对输入多项式数组执行重排操作。
 *
 * @param polOut: 输出数组,用于存储经过重排后的多项式系数。
 * @param polIn: 输入数组,包含待重排的多项式系数。
 * @param index: 重排索引数组,定义多项式系数的重排顺序。
 * @param level: 层级数量,表示要处理的多项式分块数量。
 */
void Context::PermuteNTTWithIndex(uint64_t *polOut, uint64_t *polIn, uint64_t *index, int level)
{
#pragma omp parallel for
    for (long l = 0; l < level; l++) {
        uint64_t *polIni = polIn + (l << logN);
        uint64_t *polOuti = polOut + (l << logN);
        IPermuteNTTWithIndex(polOuti, polIni, index);
    }
}

/*
 * @brief :使用指定的重排索引数组,对输入多项式数组进行重新排列。
 *
 * @param polOut: 输出数组,用于存储重排后的多项式系数。
 * @param polIn: 输入数组,包含待重排的多项式系数。
 * @param index: 重排索引数组,指定多项式系数的重排顺序。
 */
void Context::IPermuteNTTWithIndex(uint64_t *polOut, uint64_t *polIn, uint64_t *index)
{
    for (int i = 0; i < N; i++) {
        polOut[i] = polIn[index[i]];
    }
}

/*
 * @brief :对整数数组进行编码,将数据映射到多项式系数表示中。
 *
 * @param a :输出的多项式数组,每个槽的整数值将被编码为整数并存储在该数组中。
 * @param v :输入的整数数组,包含待编码的数据,长度为 slots。
 * @param slots :表示数据的槽数,决定输入数组的大小。
 * @param l :编码过程中使用的模数个数,用于多项式系数的映射。
 */
void Context::Encode(uint64_t *a, uint64_t *v, long slots, long l)
{
    for (long i = 0; i < slots; i++) {
        b[indexMatrix[i]] = v[i];
    }
    for (uint64_t j = 0; j < N; j++) {
        ModBarrett(b[j], b[j], tVec[0], twotVec[0], trVec_high[0], trVec_low[0]);
    }
    for (uint64_t i = slots; i < N; i++) {
        b[indexMatrix[i]] = 0;
    }
    TiINTTAndEqual(b, 0);
#pragma omp parallel for
    for (uint64_t i = 0; i < l; i++) {
        std::copy(b, b + N, a + i * N);
        uint64_t *ai = a + i * N;
        QiNTTAndEqual(ai, i);
        QiMulConstAndEqual(ai, t0InvModQ[i] % qVec[i], i);
    }

}

/*
 * @brief :对整数多项式进行解码,将其转换回整数形式的数据表示。
 *
 * @param a: 输入的多项式数组,包含经过编码和变换的数据。
 * @param v: 输出的整数数组,解码后的整数值将存储在此数组中。
 * @param slots: 数据槽的数量,表示整数数组的有效长度。
 * @param l: 多项式编码过程中使用的模数个数。
 */
void Context::Decode(uint64_t *a, uint64_t *v, long l)
{
    INTTAndEqual(a, l);
    MulConstAndEqual(a, tVec[0], l);
    if (l > 1) {
        ZZ c = qProduct[l - 1] >> 1;
#pragma omp parallel for
        for (long i = 0; i < l; i++) {
            uint64_t *ai = a + (i << logN);
            QiAddConstAndEqual(ai, c % qVec[i], i);
        }
        BaseConvertQ2T(a, l, b);
        for (uint64_t j = 0; j < N; j++) {
            SubMod(b[j], b[j], c % tVec[0], tVec[0]);
        }
    } 
    else {
        AddConstAndEqual(a, qVec[0] >> 1, l);
        for (uint64_t j = 0; j < N; j++) {
            ModBarrett(b[j], a[j], tVec[0], twotVec[0], trVec_high[0], trVec_low[0]);
        }
        for (uint64_t j = 0; j < N; j++) {
            SubMod(b[j], b[j], qVec[0] >> 1, tVec[0]);
        }
    }
    TiNTTAndEqual(b, 0);
    for (uint64_t j = 0; j < N; j++) {
        v[j] = b[indexMatrix[j]];
    }
}

/*
 * @brief :将输入多项式数组从模数 Q 分组转换到目标模数 T
 *
 * @param a: 输入数组,包含模数集合 Q 下的多项式系数。
 * @param length: 模数 Q 的层级数量。
 * @param res: 输出结果向量,存储每个分块转换后的结果数组指针。
 */
void Context::BaseConvertQ2T(uint64_t *a, long length, uint64_t *res)
{
#pragma omp parallel for
    for (auto n = 0; n < N; n++) {
        uint64_t ra[BASECONVERT];
        float vv = 0;
        for (auto j = 0; j < length; j++) {
            uint64_t *aj = a + (j << logN);
            MulModShoup(ra[j], aj[n], qHatInvModq[length - 1][j], qHatInvModqShoup[length - 1][j], qVec[j]);
            vv += ra[j] / static_cast<float>(qVec[j]);
        }
        for (auto j = 0; j < T; j++) {
            __uint128_t sum = 0;
            for (auto m = 0; m < length; m++) {
                __uint128_t tt = ra[m];
                sum += tt * qHatModt[length - 1][m][j];
            }
            unsigned __int128 QMul = static_cast<unsigned __int128>(vv) * QModt[length - 1][j];
            ModBarrett(res[(j << logN) + n], sum - QMul, tVec[j], twotVec[j], trVec_high[j], trVec_low[j]);
        }
    }
}

/* *
 * @brief 执行密文的小工具乘积运算，包括基础转换和模数降级
 * 当levelP>1时执行多层级计算，否则执行单层级计算和位分解
 * 最后通过并行方式进行两个基础转换操作，完成模数转换
 *
 * @param levelQ 模数Q的层级数量
 * @param levelP 模数P的层级数量
 * @param cx 输入密文数组
 * @param gadgetCiphertext 小工具密文密钥
 * @param ct 输出密文结果
 */
#include "timer_utility.h"
void Context::GadgetProduct(int levelQ, int levelP, uint64_t *cx, SwitchKey &gadgetCiphertext, uint64_t *ax, uint64_t *bx)
{
    if (levelP > 1) {
        GadgetProductNoModDown(levelQ, levelP, cx, gadgetCiphertext, tmpAx, tmpBx);
    } else {
        GadgetProductSinglePAndBitDecompNoModDown(levelQ, levelP, cx, gadgetCiphertext, tmpAx, tmpBx);
    }
    BaseConvertQlPk2Ql(ax, tmpAx, levelQ, levelP);
    BaseConvertQlPk2Ql(bx, tmpBx, levelQ, levelP);
}

/* *
 * @brief 执行不带模数降级的小工具乘积运算，包含RNS分解和多项式乘法
 * 首先对输入密文执行逆NTT变换，然后在RNS基础上分解处理，
 * 通过循环累加方式计算最终的乘积结果
 *
 * @param levelQ Q基数的层级数量
 * @param levelP P基数的层级数量
 * @param cx 输入密文数组
 * @param gadgetCiphertext 用于小工具乘积的密钥
 * @param tmpAx 存储中间计算结果的第一部分
 * @param tmpBx 存储中间计算结果的第二部分
 */
void Context::GadgetProductNoModDown(int levelQ, int levelP, uint64_t *cx, SwitchKey &gadgetCiphertext, uint64_t *tmpAx,
    uint64_t *tmpBx)
{
    INTT(cxInvNTT, cx, levelQ);
    int beta = DecompRNS(levelQ, levelP);
    for (int i = 0; i < beta; i++) {
        DecomposeSingleNTT(levelQ, levelP, i, cx, cxInvNTT, c2QP, c2QP + levelQ * N);
        if (i == 0) {
            Mul(tmpAx, gadgetCiphertext.rax[i][0], c2QP, levelQ);
            for (int j = 0; j < levelP; j++) {
                PiMul(tmpAx + (levelQ + j) * N, gadgetCiphertext.rax[i][0] + (L + j) * N, c2QP + (levelQ + j) * N, j);
            }
            Mul(tmpBx, gadgetCiphertext.rbx[i][0], c2QP, levelQ);
            for (int j = 0; j < levelP; j++) {
                PiMul(tmpBx + (levelQ + j) * N, gadgetCiphertext.rbx[i][0] + (L + j) * N, c2QP + (levelQ + j) * N, j);
            }
        } else {
            Mul(polAx, gadgetCiphertext.rax[i][0], c2QP, levelQ);
            for (int j = 0; j < levelP; j++) {
                PiMul(polAx + (levelQ + j) * N, gadgetCiphertext.rax[i][0] + (L + j) * N, c2QP + (levelQ + j) * N, j);
            }
            Mul(polBx, gadgetCiphertext.rbx[i][0], c2QP, levelQ);
            for (int j = 0; j < levelP; j++) {
                PiMul(polBx + (levelQ + j) * N, gadgetCiphertext.rbx[i][0] + (L + j) * N, c2QP + (levelQ + j) * N, j);
            }
            AddAndEqual(tmpAx, polAx, levelQ, levelP);
            AddAndEqual(tmpBx, polBx, levelQ, levelP);
        }
    }
}

/* *
 * @brief 执行单个RNS分量的NTT变换处理
 * 对密文进行分解和分割后，分别在Q基和P基上执行NTT变换，
 * 处理特定范围内的系数并进行变换
 *
 * @param levelQ Q基数的层级数量
 * @param levelP P基数的层级数量
 * @param decompRNS RNS分解的索引
 * @param c2NTT 输入的NTT域系数
 * @param c2InvNTT 输入的逆NTT域系数
 * @param c2QiQ Q基上的输出结果
 * @param c2QiP P基上的输出结果
 */
void Context::DecomposeSingleNTT(int levelQ, int levelP, int decompRNS, uint64_t *c2NTT, uint64_t *c2InvNTT,
    uint64_t *c2QiQ, uint64_t *c2QiP)
{
    DecomposeAndSplit(levelQ, levelP, decompRNS, c2InvNTT, c2QiQ, c2QiP);
    int p0idxst = decompRNS * levelP;
    int p0idxed = p0idxst + levelP;
#pragma omp parallel for
    for (int x = 0; x < levelQ; x++) {
        if (p0idxst <= x && x < p0idxed) {
            std::copy(c2NTT + x * N, c2NTT + (x + 1) * N, c2QiQ + x * N);
        } else {
            QiNTTAndEqual(c2QiQ + N * x, x);
        }
    }
#pragma omp parallel for
    for (int i = 0; i < levelP; i++) {
        uint64_t *c2QiPi = c2QiP + i * N;
        PiNTTAndEqual(c2QiPi, i);
    }
}

/* *
 * @brief 执行RNS分解和分割操作
 * 根据Q基和P基的层级数量计算分解层级，并处理系数转换。
 * 当分解层级为-1时直接复制，否则执行基于QBig的模转换
 *
 * @param levelQ Q基数的层级数量
 * @param levelP P基数的层级数量
 * @param decompRNS RNS分解的索引
 * @param p0Q 输入的Q基系数数组
 * @param p1Q 输出的Q基转换结果
 * @param p1P 输出的P基转换结果
 */
void Context::DecomposeAndSplit(int levelQ, int levelP, int decompRNS, uint64_t *p0Q, uint64_t *p1Q, uint64_t *p1P)
{
    int lvlQStart = decompRNS * levelP;
    int decompLvl = 0;

    if (levelQ > levelP * (decompRNS + 1)) {
        decompLvl = levelP - 2;
    } else {
        if (levelP == 0) {
            return;
        }
        decompLvl = ((levelQ - 1) % levelP) - 1;
    }
    if (decompLvl == -1) {
        for (int j = 0; j < levelQ; j++) {
            std::copy(p0Q + lvlQStart * N, p0Q + (lvlQStart + 1) * N, p1Q + j * N);
        }
        for (int j = 0; j < levelP; j++) {
            std::copy(p0Q + lvlQStart * N, p0Q + (lvlQStart + 1) * N, p1P + j * N);
        }
    } else {
        int p0idxst = decompRNS * levelP;
        int p0idxed = p0idxst + levelP;
        if (p0idxed > levelQ) {
            p0idxed = levelQ;
        }

        ZZ QBig = conv<ZZ>(1);
        for (int i = p0idxst; i < p0idxed; i++) {
            QBig *= qVec[i];
        }
        ZZ QHalf = QBig >> 1;
        uint64_t *QHalfModqi = new uint64_t[p0idxed - p0idxst];
        for (int i = 0, j = p0idxst; j < p0idxed; i++, j++) {
            QHalfModqi[i] = QHalf % qVec[j];
        }

        SlicesBaseConvertQl2QlP(p0Q, QHalfModqi, p1Q, p1P, levelQ, levelP, decompLvl, p0idxed - p0idxst, decompRNS,
            p0idxst, QHalf);
    }
}

/* *
 * @brief 将输入多项式数组从模数 Q 转换到目标模数 QP
 *
 * @param a
 * 输入的多项式数组，该数组来自于整个模数Q上的数组按模数P的层数进行拆分后得来
 * @param QHalfModqi 模数Q的一半在每一层上的余数
 * @param p1Q 转换后的多项式在模数Q上的部分
 * @param p1P 转换后的多项式在模数P上的部分
 * @param levelQ 模数Q的层级
 * @param levelP 模数P的层级
 * @param decompLvl 等于length - 2
 * @param length 当前待转换的多项式在模Q上的层数
 * @param slices 当前待转换的多项式的分组索引
 * @param offset 当前待转换的多项式的在原模Q数组上的偏移量
 * @param QHalf 模数Q的一半
 */
void Context::SlicesBaseConvertQl2QlP(uint64_t *a, uint64_t *QHalfModqi, uint64_t *p1Q, uint64_t *p1P, int levelQ,
    int levelP, int decompLvl, long length, long slices, long offset, ZZ QHalf)
{
#pragma omp parallel for
    for (auto j = 0; j < length; j++) {
        uint64_t *aj = a + ((j + offset) << logN);
        uint64_t *tmpj = tmpQl2QlP + (j << logN);
        for (auto n = 0; n < N; n++) {
            MulModBarrett(tmpj[n], aj[n] + QHalfModqi[j], qInvModqWithQP[levelP - 2][slices][decompLvl][j],
                qVec[offset + j], twoqVec[offset + j], qrVec_high[offset + j], qrVec_low[offset + j]);
        }
    }
#pragma omp parallel for
    for (auto j = 0; j < levelQ; j++) {
        if (j >= offset && j < offset + length) {
            continue;
        }
        for (auto n = 0; n < N; n++) {
            uint64_t tt = tmpQl2QlP[n];
            unsigned __int128 sum =
                static_cast<unsigned __int128>(tt) * qInvModqpWithQP[levelP - 2][slices][decompLvl][j][0];
            double v = tt / (qVec[offset] * 1.0);
            for (auto m = 1; m < length; m++) {
                tt = tmpQl2QlP[n + (m << logN)];
                sum += static_cast<unsigned __int128>(tt) * qInvModqpWithQP[levelP - 2][slices][decompLvl][j][m];
                v += tt / (qVec[offset + m] * 1.0);
            }
            unsigned __int128 QMul = static_cast<unsigned __int128>(v) * QModqpWithQP[levelP - 2][slices][decompLvl][j];
            uint64_t QHalfMod = QHalf % qVec[j];
            while (sum - QMul < 0) {
                sum += qVec[j];
            }
            sum -= QMul;
            if (sum - QHalfMod < 0) {
                sum += qVec[j];
            }
            sum -= QHalfMod;
            ModBarrett(p1Q[(j << logN) + n], sum, qVec[j], twoqVec[j], qrVec_high[j], qrVec_low[j]);
        }
    }
    for (auto j = L; j < L + levelP; j++) {
        for (auto n = 0; n < N; n++) {
            uint64_t tt = tmpQl2QlP[n];
            unsigned __int128 sum =
                static_cast<unsigned __int128>(tt) * qInvModqpWithQP[levelP - 2][slices][decompLvl][j][0];
            double v = tt / (qVec[offset] * 1.0);
            for (auto m = 1; m < length; m++) {
                tt = tmpQl2QlP[n + (m << logN)];
                sum += static_cast<unsigned __int128>(tt) * qInvModqpWithQP[levelP - 2][slices][decompLvl][j][m];
                v += tt / (qVec[offset + m] * 1.0);
            }
            unsigned __int128 QMul = static_cast<unsigned __int128>(v) * QModqpWithQP[levelP - 2][slices][decompLvl][j];
            uint64_t QHalfMod = QHalf % pVec[j - L];
            while (sum - QMul < 0) {
                sum += pVec[j - L];
            }
            sum -= QMul;
            if (sum - QHalfMod < 0) {
                sum += pVec[j - L];
            }
            sum -= QHalfMod;
            ModBarrett(p1P[((j - L) << logN) + n], sum, pVec[j - L], twopVec[j - L], prVec_high[j - L],
                prVec_low[j - L]);
        }
    }
}

/* *
 * @brief 执行单P基数下的小工具乘积和位分解运算，不包含模数降级
 * 首先执行逆NTT变换，然后在Q基和P基上分别进行NTT域乘法运算，
 * 并通过循环累加方式合并结果
 *
 * @param levelQ Q基数的层级数量
 * @param levelP P基数的层级数量，在此函数中为单层
 * @param cx 输入密文数组
 * @param gadgetCt 用于小工具乘积的密钥
 * @param tmpAx 存储结果的第一部分
 * @param tmpBx 存储结果的第二部分
 */
void Context::GadgetProductSinglePAndBitDecompNoModDown(int levelQ, int levelP, uint64_t *cx, SwitchKey &gadgetCt,
    uint64_t *tmpAx, uint64_t *tmpBx)
{
    INTT(cxInvNTT, cx, levelQ);
    memset(tmpAx, 0, sizeof(uint64_t) * (levelQ + levelP) * N);
    memset(tmpBx, 0, sizeof(uint64_t) * (levelQ + levelP) * N);
    #pragma omp parallel for
    for(int i = 0; i < levelQ; i++) {
        for(int j = 0;j < levelQ; j++) {
            std::copy(cxInvNTT + i * N, cxInvNTT + (i + 1) * N, cwQSet[i] + j * N);
        }
        std::copy(cxInvNTT + i * N, cxInvNTT + (i + 1) * N, cwPSet[i]);
        NTTAndEqual(cwQSet[i], levelQ);
        PiNTTAndEqual(cwPSet[i], 0);
    }
    for (int i = 0; i < levelQ; i++) {    
#pragma omp parallel for collapse(2)
        for (int u = 0; u < levelQ; u++) {
            for (long n = 0; n < N; n++) {
                MulModBarrett((polAx + u * N)[n], (gadgetCt.rax[i][0] + u * N)[n], (cwQSet[i] + u * N)[n], 
                                    qVec[u], twoqVec[u], qrVec_high[u], qrVec_low[u]);
                MulModBarrett((polBx + u * N)[n], (gadgetCt.rbx[i][0] + u * N)[n], (cwQSet[i] + u * N)[n], 
                                    qVec[u], twoqVec[u], qrVec_high[u], qrVec_low[u]);
                AddMod((tmpAx + u * N)[n], (tmpAx + u * N)[n], (polAx + u * N)[n], qVec[u]);
                AddMod((tmpBx + u * N)[n], (tmpBx + u * N)[n], (polBx + u * N)[n], qVec[u]);
            }
        }
        PiMul(polAx, gadgetCt.rax[i][0] + L * N, cwPSet[i], 0);
        PiMul(polBx, gadgetCt.rbx[i][0] + L * N, cwPSet[i], 0);
        PiAddAndEqual(tmpAx + levelQ * N, polAx, 0);
        PiAddAndEqual(tmpBx + levelQ * N, polBx, 0);
    }
}


/* *
 * @brief 将输入多项式数组从模数 QP 转换到目标模数 Q
 *
 * @param res 存储转换后的多项式数组
 * @param a 输入的待转换的多项式数组
 * @param levelQ 模数Q的层级
 * @param levelP 模数P的层级
 */
void Context::BaseConvertQlPk2Ql(uint64_t *res, uint64_t *a, long levelQ, long levelP)
{
    for (int i = 0; i < levelP; i++) {
        PiINTTAndEqual(a + (levelQ + i) * N, i);
    }
    ZZ pHalf = pModulusAtLevel[levelP - 1] >> 1;
    for (int i = 0; i < levelP; i++) {
        PiAddConstAndEqual(a + (levelQ + i) * N, pHalf % pVec[i], i);
    }
    for (int i = 0; i < levelP; i++) {
        uint64_t *ai = a + ((i + levelQ) << logN);
        uint64_t *tmpi = tmpP2Ql + (i << logN);
    #pragma omp parallel for
        for (int n = 0; n < N; n += 2) {
            uint64x2_t ai_vec = vld1q_u64(ai + n);
            uint64x2_t tmpi_vec;
            uint64x2_t pSliceHatInvModp_vec = vdupq_n_u64(pSliceHatInvModp[levelP - 1][i]);
            uint64x2_t pVec_vec = vdupq_n_u64(pVec[i]);
            MulModShoup_NEON(&tmpi_vec, ai_vec, pSliceHatInvModp_vec,
                vdupq_n_u64(x_Shoup(pSliceHatInvModp[levelP - 1][i], pVec[i])), pVec_vec);
            vst1q_u64(tmpi + n, tmpi_vec);
        }
    }
#pragma omp parallel for collapse(2)
    for (int i = 0; i < levelQ; i++) {
        for (int n = 0; n < N; n += 2) {
            int offset = i << logN;
            uint64x2_t tmp_vec = vld1q_u64(tmpP2Ql + n);
            __uint128_t sum_low =
                static_cast<__uint128_t>(vgetq_lane_u64(tmp_vec, 0)) * pSliceHatModq[levelP - 1][0][i];
            __uint128_t sum_high =
                static_cast<__uint128_t>(vgetq_lane_u64(tmp_vec, 1)) * pSliceHatModq[levelP - 1][0][i];
            double v_low = vgetq_lane_u64(tmp_vec, 0) / (pVec[0] * 1.0);
            double v_high = vgetq_lane_u64(tmp_vec, 1) / (pVec[0] * 1.0);
            for (int m = 1; m < levelP; m++) {
                tmp_vec = vld1q_u64(tmpP2Ql + n + (m << logN));
                sum_low += static_cast<__uint128_t>(vgetq_lane_u64(tmp_vec, 0)) * pSliceHatModq[levelP - 1][m][i];
                sum_high += static_cast<__uint128_t>(vgetq_lane_u64(tmp_vec, 1)) * pSliceHatModq[levelP - 1][m][i];
                v_low += vgetq_lane_u64(tmp_vec, 0) / (pVec[m] * 1.0);
                v_high += vgetq_lane_u64(tmp_vec, 1) / (pVec[m] * 1.0);
            }
            uint64_t PMul_low = static_cast<uint64_t>(v_low) * PSliceModq[levelP - 1][i];
            uint64_t PMul_high = static_cast<uint64_t>(v_high) * PSliceModq[levelP - 1][i];
            uint64_t result_low, result_high;
            ModBarrett(result_low, sum_low - PMul_low, qVec[i], twoqVec[i], qrVec_high[i], qrVec_low[i]);
            ModBarrett(result_high, sum_high - PMul_high, qVec[i], twoqVec[i], qrVec_high[i], qrVec_low[i]);
            res[offset + n] = result_low;
            res[offset + n + 1] = result_high;
        }
    }
#pragma omp parallel for
    for (int i = 0; i < levelQ; i++) {
        QiSubConstAndEqual(res + i * N, pHalf % qVec[i], i);
    }
    NTTAndEqual(res, levelQ);
#pragma omp parallel for collapse(2)
    for (int i = 0; i < levelQ; i++) {
        for (int n = 0; n < N; n += 2) {
            int offset = i << logN;
            uint64x2_t a_vec = vld1q_u64(a + offset + n);
            uint64x2_t res_vec = vld1q_u64(res + offset + n);
            uint64x2_t twoq_vec = vdupq_n_u64(qVec[i] << 1);
            uint64x2_t mid_vec = vsubq_u64(vaddq_u64(twoq_vec, vsubq_u64(res_vec, a_vec)), vdupq_n_u64(0));
            uint64x2_t result_vec;
            MulModShoup_NEON(&result_vec, mid_vec, vdupq_n_u64(qVec[i] - PSliceInvModq[levelP - 1][i]),
                vdupq_n_u64(x_Shoup(qVec[i] - PSliceInvModq[levelP - 1][i], qVec[i])), vdupq_n_u64(qVec[i]));
            vst1q_u64(res + offset + n, result_vec);
        }
    }
}

/* *
 * @brief 扩展小范数系数到P基并进行中心化处理
 * 根据输入系数与QHalf的比较确定符号，
 * 将结果扩展到levelP个P基中
 *
 * @param res 输出的扩展结果数组
 * @param a 输入的小范数系数数组
 * @param levelP P基数的层级数量
 */
void Context::ExtendBasisSmallNormAndCenter(uint64_t *res, uint64_t *a, int levelP)
{
    uint64_t QHalf = qVec[0] >> 1;
    for (int j = 0; j < N; j++) {
        uint64_t coeff = a[j];
        int sign = 1;
        if (coeff > QHalf) {
            coeff = qVec[0] - coeff;
            sign = 0;
        }
        for (int i = 0; i < levelP; i++) {
            uint64_t pi = pVec[i];
            res[i * N + j] = (coeff * sign) | (pi - coeff) * (sign ^ 1);
        }
    }
}